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A vortex model is used to develop a scheme to compute the 
pressure field around a marine propeller operating in a wa e w ic 
is nonuniform in both the radial and circumferential directiono. 

This model takes into account the effects of ° distri- 

outline of the blades and the form of the unsteady loading distri 
bution, both radially and over the chord. The model can be used 
for lightly as well as "moderately" loaded propellers. Following 
the theoretical development a computer program is written in 
fortran IV to be used with the IBM System/360-65 digital compute^ . 
S ^habUity of the program Is tested with sample calculations. 
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Chapter I 
INTRODUCTION 

The propeller induced vibratory forces acting on a ship have 
been for a long time an important subject of experimental and theo- 
retical research. Although some important conclusions have been 
reached based on reliable experimental results [17], it is always 
hard, if not impossible, to support such conclusions with theoretical 
calculations. One reason for this defficiency is that for a long time 
no analytical method was available for evaluating the propeller loading 
as a three-dimensional problem. 

The initial numerical results obtained for the propeller-induced 
pressure field were based on a simple lifting-line representation [lO] 
which was found to predict pressures far below those obtained from 
measurement. This theory was later refined by including the effect of 
the blade thickness and this brought the numerical results closer to 
the experimental data [9, 10, 11]. 

The effects on the pressure field due to nonuniform inflow were 
studied by Tsakonas ^ a_l. [8], They treated the problem as a Dirichlet 
problem of potential theory and related the pressure field near the 
propeller to the trust- and torque-producing forces. These forces 
were computed using the gust theory for airfoil^ and assuming that 
the flow around each blade was two-dimensional. Two correction factors 
were added to compensate for cascade and finite aspect ratio effects. 
Even though the calculations made using this method cannot be con- 
sidered very accurate, they can give an idea of the relative importance 
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of the unsteady effects arising from the nonuniformity of the ship^s 
wake. This can be seen in Fig. taken from [8]. It is apparent 
that very close to the plane of the propeller the increase in pressure 
due to nonuniformity is at most 20 percent of the total pressure but 
beyond x/D =0.4 the percentage change becomes much greater. The 
reason is that at such distances all the pressures are small^ but the 
slow rate of decay of the components arising from nonuniformity makes 
them more important. It must be noted^ however^ that the pressure due 
to thickness is not included in Fig. 1 so that the relative effect 
close to the propeller arising from nonuniform inflow is of the order 
of 10 percent for the particular wake used in these calculations. 

A truly three-dimensional theory was recently used by Kerwin 
[18] to compute the field point velocities induced by a marine pro- 
peller operating in a steady flow. He makes use of the vortex theory 
to represent the blade loading and derives a scheme to compute the per- 
turbation velocities on specified field points located on the blades. 

In the present work we will extend the program developed by 
Professor Kerwin to include field points located at any point in the 
fluid around the propeller and to take into account the unsteady effects 
due to a nonuniform inflow such as is found in a ship^s wake. In 
addition^ a scheme will be developed to compute the velocity potential 
at the field points. The pressure field is then computed by sub- 
stituting the perturbation velocities and time variation of the 
velocity potential into Bernoulli’s equation. 
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FIG. 1 ~ PRESSURE DISTRIBUTION FOR THE FIRST BLADE HARMONIC OF A 4 BLADE PROPPELER 
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The usual assumptions of an incompressible^ frictionless and 
irrotational flow are held. It is also assumed that all the disturb- 
ances are small, which requires that the blades be thin and propeller 
loading small. The effects of thickness form and outline of the 
blades and the form of the loading distribution, both radially and 
over the chord are taken into account by this theory. 

When this work was already under way, the Davidson Laboratory 
published a report by Jacobs al. [19] which also treates the un- 
steady, propeller induced, pressure field using a lifting surface 
theory. However, their treatment is based on the solution of a Dirichlet 
problem for the pressure field. It seemed of interest to compare the 
results obtained with the two different approaches but unfortunately 
were were not able to obtain from D.L, the loading distribution for 
the propellers used in their report. 
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Chapter II 



FORMULATION OF THE PROBLEM 



It is desired to find the pressure at any point in the fluid 
around a ship*s propeller. Assuming that the flow is incompressible, 
frictionless and irrotational, the problem becomes one of potential 



potential and velocities induced by the propeller and substituting 
into the Bernoulli's equation to solve for the pressure. 

^ • Geometry of the Propeller 

The following geometrical considerations are adopted from [!]• 

Two orthogonal coordinate systems -- a Cartesian (x,y,z) and a 
cylindrical (x,r,6) -- are associated with the propeller as represented 
in Fig. 2. Both systems are fixed with the propeller and the first 
blade is made to coincide with the y axis. All length dimensions are 
nondimensionalized with respect to the radius of the propeller, R, 

The equations relating the two coordinate systems are: 



In order to relate corresponding points on each of the K blades, 
we define as the 6 coordinate of the point at the tip of the k'th 
blade. For symmetrically arranged blades it is 



theory which in principle can be solved by finding the velocity 




( 1 ) 



= 2TT(k- 1) 
k K 




( 2 ) 
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FIG. 3 - PROJECTED BLADE OUTLINE 
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This is illustrated in Fig. 3 which also shows the projected blade 
outline, defined by the angular coordinates 6 ^ (r) and 0 (r) of the 
leading and trailing edges, respectively. 

The geometry of the flow past a blade is schematically re- 
presented in Fig. 4. It is seen that the oncoming flow V^(r) forms 
an angle P with the yz plane, and we define 

X(r) = r tanp(r) (3) 

as the ’’advance coefficient”. The streamlines are helices with 
pitch 

P(r) = 2v X(r) • (4) 

It is useful to define an orthogonal curvilinear coordinate 
system (s,r,n) associated with the helical surfaces swept out by 
the undisturbed flow past the radial lines {x = 0, 0 = This 

system, shown in Fig. 2^ has unit vectors; 

u =i sinP - sin(0 +6, ) cosP +k cos (0 +6, ) cosP 

S K K 

= cos(6 +5j^) +k sinC© (5) 

u =i cosp + sin(0 +6 ) sinP - ic cos (0 +6, ) sinp 
n K. rC 

It is assumed that the blades lie on those helical surfaces. 
This assumption is in accordance with the linear theory for lightly 
loaded propellers, however, it does not impose any limitation to the 
possibility of including nonlinear refinements by suitably modifying 
the advance angle [*1]. 
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FIG. 4 - APROACH FLOW \^LOCITY DIAGRAM 
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2 . Expression for the Pressure 

The Bernoulli's equation for unsteady flow is 

2 2 

P = ®f ^ 2 - V ) 

where 

p = hydrodynamic pressure 
= fluid density 
$ = velocity potential 

= speed of free stream or speed of advance 
V = resultant velocity at a point 
t = time variable. 

Changing to a frame of reference moving with the propeller 
blades and after linearization the above equation becomes [11] 



P ■ ^ PfVA Pf''' "t 



( 6 ) 



where 



u ^u = axial and tangential disturbance velocities 
a t 

Q = angular velocity of the propeller. 



3 . Induced Velocities 

In order to find the velocities induced by the propeller^ we 
make use of the vortex representation described^ for instance^ in 
[1] and [2]. The following is an adaptation to the unsteady case 
of the relations derived in [1]. 



i 



11 - 



The inflow velocity and verticity strength are expressed as 
complex Fourier series of the angular position of the rotating co- 
ordinate system with respect to a set of axes fixed in the ship^ 

v^(p,cpt) = V^(p) 

r(p^cp,t) = r(p, cp) -e 

where 

V = relative flow velocity past the blade section 
o 

X = vorticity strength^ in general 
N = harmonic number 

Q = angular velocity of the propeller 



(i) Velocity Induced by Source Distribution 

The velocity induced at a point (x^r^0) by sources distributed 
over K blades is 






K 

= E 

k=l 



^iN(^2t-5k) 

kir 



.1 0n 



n 



V„(p) 



. [x-Mp)-9i_[p^+A^(p)f ^ jp (7) 



K 



,e ,t) =Tj 

k=l 






kv 




A ®T 



Vo(p) 



iNp OT 
• e ^ • ‘c — • 
os 



P sin(cp+ - e) • [p^+ A^(p) ] 



2 . -, 2 , , 



dcp dp (8) 
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where 



D = [(x -Mp) *cp) + r + p - Zrp cos(cp+5^- e; 




X = blade thickness 



The superscript s denotes a source velocity. 

(ii) Velocity Induced by Vortex Distribution 
We separate the components of the induced velocities due to 
the bound vorticity in the blade (superscript b) ^ radial vorticity 
in the wave (sup. r)^ streamwise vorticity in the blade (sup. f) 
and streamwise vorticity in the wake (sup. w) . 




h 



r[p^+A^(p)]^ sin(cp+ e) 



dcp dp 



(9) 





dcp dp (10) 




r[p^+A^(p)]^ sin(cp+ 6j^ - 0) 



( 11 ) 
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(x,r,e,t) = Tj 



k-l ri, 01 



(x- A(p)(p)[p^+ A^(cp)]^ cos(tp+ 6j^ - 0) 



dcp dp (12) 



iN(fit-5 ) 1 e 

^ iV 0 I I -L 

“a fe,r, 0 ,l:) -p 5 ^ 

Ph H 



fl' 






[p - pr cos(cp+ 5, - 0) ] 



dq) dp 



(13) 



K 



u^^^(x,r,e,t) =S j— 

k=l 



iN(fit-5, ) -1 



n 



f.(p,cp) 



{ 



?\(p)[r- p cos(cp+6, - 0)] 



p[x-X(p)cp] sin(cp+ 5, - 0)^ 

+ 3 ^ \ dtp dp 



(14) 



u<“>(x,r,0,t) =S 

“ k-1 



K .1 



hiT 



i Ivp)-p'^'^. 






[p - pr cos(cp+5. - 0)] 



dcp dp 



(15) 



“t“ (*.P.9.P) -p 55 p ^ w 

Ph ei 






{ 



^(p)C’^" P cos(i-p+5, - 0) p[x-'A(p)cp]sin(cp+ 6, - 0) 



D' 



•j- dp dp 



( 16 ) 
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where 



T^(p, cp) = complex amplitude of the bound vorticity 

“ complex amplitude of the radical vorticity 



in the wake 



T^Cp^Cp) = complex amplitude of the streamwise vorticity 
in the blade 



“ complex amplitude of the streamwise vorticity 



in the wake. 



3. Velocity Potential 

As before we separate the effects due to source and vortex 
distributions . 

(i) Velocity Potential due to Source Distribution 
The velocity potential at a point due to a source of unit 
strength located at a distance^ d^ is given by [16] 



Multiplying by the source strength integrating over the first 
blade and summing over the K blades gives: 






( 17 ) 



( 

I 



i 

1 



- 15 - 



(ii) Velocity Potential due to Vortex Distribution 
Instead of trying to find an exact expression for the velocity 
potential^ which would be rather cumbersome^ we make use of the fact 
that the unsteady effects become important only at a distance from 
the propeller and that a reentrant vortex is^ in a sense^ equivalent 
to a uniform distribution of double sources over any surface bounded 
by it. It is shown in Appendix A that if we substitute a reentrant 
vortex by an equivalent concentrated doublet^ the error in the value 
of the velocity potential decreases rapidly with distance. 

In this way we can compute the velocity potential by summing 
a set of discrete values given by expressions of the form (see Ap- 
pendix A) 



S cos Y 

47T d 



where 



S = strength of the concentrated doublet 
Y = angle between the doublet axis and the line from 
the doublet to the field point. 

The value of S will be found later and cos Y is easily computed 

from 



cos Y = 



n* D 

Inl-D 



where 



n = vector direction of the double axis 
D = [x - A(p)cp]i+ [r COS0 - p cos(cp +5j^) ] j + 
+ [r sin0 - p sin(cp+ 6j^) ]k 



( 



I 



I 

I 



- 16 - 



Using (5) to obtain n^ it follows that 

[x - >\(p)cp]cosP(p) + r 



COS T = 



D 



sinp(p)cos(cp+ 5^ “ 0) 



k . The Loading Modes 

For the purposes of this work it is assumed that the distri- 
bution of the loading on the propeller blades has already been de- 
termined by any of the known methods (examples [2]^ [3], [*7], [12] 
and [15]) which solve the unsteady lifting surface problem of a 
propeller operating in a wake. The usual procedure in these methods 
is to assume that the loading may be expressed in separable form as 
the product of a function of the spanwise coordinates and a function 
of the chordwise coordinates. In selecting these functions^ one 
makes use of the known results of the theory for incompressible flow 
in two dimensions and of the lifting line theory for the higher 
aspect ratio wings. This leads to functions of the type 



where 



Ap = S F^(s) X G^(r) 
_ n n 

n=0 



Ap = loading per unit area 
F (s) = chordwise function 



G^(r) = spanwise function. 



Two dimensional airfoil theory suggests that we choose F of 



n 



the form, [5] 



F (1^ = r - — 



+ sin[n(cos'^ I) ] 



(19) 
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where 

2 (s - s^e) 



I - chord length 



Introducing a new coordinate Q defined by 



s = 2 — ^ ‘ 2 0 < 9 < ir 



F becomes 
n 



/S 0 /s 

F (0) = cotan -r + cos (n7T) sin(n0) 



The spanwise function is either determined directly with 

the aid of a specially selected interpolation function [5] or ex- 



pressed in the form 



G 

n 



(r) = S A H (r) 
nm m 



m 

and then determined by the solution of the integral equation resulting 
from the formulation of the lifting surface problem. In order to ob- 
tain a numerical solution G^(r) is usually approximated by a stepwise 
distribution. 

The loading then will be of the form 



n °° 

Ap(r^^= L (r) cot + Z) L (r) sin(n0) (20) 

o Z ' n 

n=l 



where L and L are stepwise functions of the radius. This is called 
on ^ 

the Birnbaum distribution. 



When the loading is expressed in terms of the bound vorticity^ 



it is necessary to add an extra term to the Birnbaum modes in order 
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to satisfy the Kutta condition at the trailing edge [13]. In this 
case it becomes 

00 

r^(r,0) = A^(r) cot | + Z/ A^(r) sin(ne) + f(r^9) (21) 

n=l 

where the Kutta condition requires that 



f(r,7T) =r^(r,TT) = -ik^F (22) 

where 

_ N n £(r) 

n ~ 2V (r) 
o 

r = chordwise integration of 
According to Brown [3] a convenient choice for f(r^0) is 

/N 

.1 T- 1 - COS0 

f(r,0) = X 2 

After substituting into (21) and integrating over the chord^ it fol- 
lows that 



f(r,e) = -TT X i-jTTiT X i X 



(23) 



It must be noted that there is no correspondence (except for 
the zeroth harmonic) between the coefficients of (20) and (21) . 
However^ when the loading is given in terms of pressure distribution 
one can always find the corresponding vorticity distribution from 
the solution of the following Volterra equation of the second kind 



[14] 



^ ^ i rj^(r,s,t)ds + V^(r) Tb(r,s,t) 



( 24 ) 



in terms of 9 the above expression becomes 
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rjr,0) = 



Ap(r,e) 
Pf V^(r) 



r • 

\ Sll 



/\ y\ 



ik (r) \ sin0 T, (r,0) d0 
n D ^ 

0 



(25) 



5 • The Vorticity Relations 

The pulsating bound vorticity is expressed in the form 



= fu(p,q5) e 






(26) 



All the other components of the vorticity can be expressed in 
terms of as required to preserve continuity of vorticity (see [1] 
and [2] for detailed derivation) 



= f.(p,cp) e 






(27a) 



ff(p,cp) = - ^ (rt,(p,cp)-[p^ + ^^p)]^)d? 



(27b) 



r^j(p^ 9 ^k,t) = f„(p) e 



iN[nt-(cp-0^)-6j^] 



(28a) 



f - ^r(p) 



(28b) 



rj.(p,cp.k,t) = fj.(p) e 



iN[fit-(cp-0^)-6j^] 



(29a) 



where 



^ = iNr(p) 

T “* O O J 

[p +V(p)]‘ 



(29b) 



r(p) = j fb(p,cp)-[p^+ ^^(p)]" dp 



(30) 
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Notice that and (streamwise and radial vorticity in the wake) 
are expressed as the product of a complex amplitude at the trailing 
edge and a sinusoidal time and distance variation along the helical 
surfaces in the streamwise direction. This form had already been 
implied in eqs. (6)^ (7)^ (10) and (11). 
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Chapter III 
NUMERICAL SOLUTION 

The following is an outline of the computational model used 
to solve the equations derived in the previous chapter. 

The blades and wake are divided into chordwise panels and 
these panels are subdivided by radial grid lines as shown in Fig. 6. 
The scheme for dividing the span is shown below. 



dA d/2 



d 




FIG. 5 - SPAN SUBDIVISION FOR NUMERICAL SOLUTION 

The interval is first divided into MN equal spaces. The intervals 
at each end are then subdivided into half and quarter spaces as 
shown. The result is MT = MN + 4 chordwise panels. The chord lengthy 
pitchy rake^ bound circulation strength and source strength are as- 
sumed to be independent of radius within each of the panels^ with 
values matching those of the continuous functions at the mid-radius. 
The continuous trailing vortex sheet then becomes a set of MT+ 1 
concentrated trailing vortex lines at each of the panel boundaries. 
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FIG. 6 - BLADE AND WAKE LATTICE FOR NUIffiRICAL SOLUTION 
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The bound circulation in the blade and the radial circulation 
in the wake are assumed concentrated at each of the radial grid 
lines . 



The above assumptions allow us to write the integrals (7) to 
(12) as summations^ as follows: 



MT I 



(s)/ ^ iNfit ^ -iN:p(i) SB(m,i) ^ ^ ^ 

(x,r,9,t) = e 1j 1j e ^ ^ • r -2 — [x - A(m) -cp (i) ] 

m=l 1=1 



k=l D 



m 



(31) 



MT I 



t m=l i=l ^ 



K -iN5, 
Se ' 
k=l 



^m+1 



sin(cp(i) +6,-0) 






Pm 

m 



(32) 



,t) = -e 
a 



MT T 

iNfit >o ^ GB(m,i) 

•r* Li Zj j — 

m=l i=l 



, ^nri-1 

k . .X I dp 



K -iN5 

S e • sin(cp(i)+ 6j^- 0) ^ -j 

Pm ^ 



k=l 



m 



(33) 



(b)/ ^ GB(m,i) , ^ 

Uj. (x,r,0,t) = -e S S 2 ^ — • (x - A(m) -cp (i) ) 

m=l i=l 



K 



Pm+l 



•iK5 



Z) e -cos(cp(i) + 6, - 0) 
k=l ^ 



r dp 

■jp ? 



(34) 



-r-S S 

^ m=l i=I 2 



V '*1= . , , r"*' <*p 

2j e •sin(cp(i) + 6j^- 0) -J — 



k=l 



P D 
m 



( 35 ) 
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u 



(r) 



(x,r,0,t) 



e 21/ 2-/ e 
m=l i=I 

K "iN5, 

Zj e • cos(cp(i)+6. 

k=l ^ 



GR(m) 

2 



- 0 ) ,• 



• (x - A(m) -cpCi)) 




( 36 ) 



where 

SB(m^i) = amplitude of the equivalent source strength^ 

concentrated at the i’th grid line in the m^th panel. 

GB(m_,i) = amplitude of the equivalent bound vorticity, con- 
centrated at the i*th grid line in the m* th panel. 

GR(m) = amplitude of the radial vorticity in the wake^ con- 
centrated at the i^th grid line in the m^th panel. 

In the above equations the source and vortex strengths are non- 
dimensionalized with respect to 27TRV^^ and the perturbation velocities 
are non-dimensionalized with respect to V^. 

The integrals that appear in the equations above are evaluated 
between the limits of integration corresponding to the boundaries 
of each panel and can be determined analytically as follows 



Pm-hl 





2 2 

r p cos(cp(i) + 5^-0) + (x- A(m) -cp(i) ) + r 

([x- A(m)cp(i)] + r^sin^(cp(i) + 6^- 0))-D 



^mrhl 



P 



m 



dp P " cos(cp(i) + 5^- e) 

[(x - A(m)cp(i))^ + r^sin^(cp(i) + 6^-0)]*D 

The contribution from the trailing vortices -- 



^m+1 



eqs. (13)-(16)-- 



raises the problem of streamwise integrals that cannot be evaluated 
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analytically. It was found by Kerwin [18] that for the zeroth 
harmonic a simple trapezoidal integration is as efficient as more 
sophisticated quadrature formulas^ except for nearby elements. In 
order to investigate the application of this argument to other har- 
monics^ we carried out the integration for a single trailing vortex 
line located at 0.75R^ using a Gaussian quadrature and a trapezoidal 
rule. The computation of the perturbation velocities was carried 
for different points around the propeller using a third harmonic of 
the loading. We found that the results obtained with the two for- 
mulas differ in all cases by less than 2 percent. However^ due to 
the pulsating character of the wake vorticity it is necessary to keep 
the intervals of integration small compared with the wavelength of 
the perturbation. This results in a very large number of computations 
for each fortex line (about 1100 for a field point in the plane of 
the propeller). A better alternative was found by substituting the 
infinite integrals by an infinite series of finite integrals and then 
applying the Euler transformation to speed up the convergency of the 
series. Basicly the m.ethod works as follows: 

Given an infinite integral of the type 




dx 



it is transformed into the infinite series 



(k+l)7r/N 




•^=0 Tt/N 



7T/N 




k=0 '^O 
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Now^ applying the Euler Transformation [21]^ the above series 
becomes 



I 



= s 

k=0 



(-1) A 



2^+1 



a 

o 



k 



= S (-1)” 

m=0 




a 



k-m 



a 



k 




•f(t + k 7r/N)dt 



Sample computations done with this method showed good convergence 
with a fairly reduced number of terms. Also^ it was found that a 
5-point Gaussian quadrature was adequate for the computation of the 
terras of the series. 

When the above transformations are applied to eqs. (15) and 
(16)^ we obtain 




,e ,t) 



iNfit 

= e 



MT+1 



m=1 



iN0^(m) 

e 



TL(tn) 

2 



K 

s 

k=1 



-iN8, 

k 

e 



■ E 

j=o 



Ij+1 



(37) 






.7T/N 






FWA(|+ 0^+ k Tr/N)d| 



/ N -XTO MT+1 iN0„(m) , K -iN5, 

(w)/ „ iNfit T TL(m) k 

u' = e • Zj e • — e 

m=1 ^ k=1 



J=1 



I3+i 



( 38 ) 






V/N 



-iNa+0 ) 

e • FWT(|+ 0^+ k ir/N)d| 



I 



I 
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where 



TL(m) = amplitude of the m' th trailing vortex^ on the wake 

cos(cp+5j^-e)] 

FWA(cp) = r 

D 

FWT(cp) = {-p^[A(m)cos(cp+ 5j^ - e) - (x - A(m) *cp) sin(cp+ 6^^- 0) ] + 

r • A(m) )/D^ 



It is clear that eqs. (37) and (38) cannot be applied to the zeroth 
harmonic of the loading because the upper limit of the integrals 
would go to infinity. Instead^ we will use the trapezoidal rule 
which was found appropriate for this case^ as mentioned before. 
Using a trapezoidal rule^ eqs. (15) and (16) become 



(w), , iN(e^-cp(i)) TL(m,i)+ TL(m,i-l) 

-^(x,r,0,t) = e 2 j 1j e • — ^ 



m=l i=I 



K ~iN5, p -r cos(cp(i)+ 6, “ 0) 

Km K 

p • Zj e ^ 

k=1 D 



(39) 



u 



(w) 

t 



(x,r,0,t) 



iNnt iNW^-<p(i)) 

= e ' 1j 2j e 
m=1 i=I 



TL(m,i)+ TL(m,i-l) 
4 



K 

S e (-p [A(m)cos(cp(i) +6, - 0) - 

k=1 ™ 

(x - A(m)cp(i) ) • sin(cp(i) +6j^ - 6)] r A(m) )/D^ 

(40) 

The contribution from the trailing vortices on the blade ~~ eqs. (13) 
and (14) -- are computed using a trapezoidal integration^ for all the 
harmonics ^ 






I; 



J 







- 28 - 



MT+1 I 



(f), _ . iNfit ‘ ^ GT(m,i) + GT(m,i-l) 

(x,r,e,t) = e ^ Lj 1 • 

m=1 i=1 

K -iNS, p - r cos(cp(i) + S, - 0) 
— km ^ 



•p 2-1 e 
m , , 

k=1 



( 41 ) 



MT+1 T 

‘ m =1 1-1 * 



K -1N6, 

•S e * 

k=1 



{-p [A(m)cos(cp(i)+ 5, - 6) - 
m K. 



(x - 'X(m)cp(i) • sin(cp(i) + - 0) ] +r A(m) }/D (42) 

where 

GT(m^i) = Amplitude of the m* th trailing vortex located 
on the blade ^ between the i^th and the i*th+1 
grid lines. 

The triple summations corresponding to the velocity potential 
due to the source and vortex distributions -- eqs. (17) and (18)-- are 






= -e 



iNfit 



E e 
k=1 



MX 


I 


E 


E e 


m=1 


i=1 


iN8, 


pP. 


k 


. 




J 

^m 



-iN:p(i) SB(m,i) 
2 



mfl 



D 



(43) 



/UN -MO. MX I K -iN8, 

ji.(^)/ ^ iNfit T 7 \ GP(in,i) k 

(x,r,0,t) = -e 2^ 1j — ’ Zy e 

m=1 i=1 k=1 

(x - A(m) *cp (i) ) cos p(m)+r sin P(m) cos(cp(i) + 5, - 0) 



( 44 ) 
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(b) 



MT M 



w 



/ ^ \ iN^t ^ V*' GL(m,i) 

(x,r,e,t) = -e 1j 1j 7^ — 2^ 



K -iN5, 



m=1 i=1 



k=1 



(x - A(m) *cp(i))cos P(m) + r sin P(m) cos (cp(i) + 6, - 0) 



(45) 



where 

= velocity potential due to the vorticity in the blade. 
= velocity potential due to the vorticity in the wake. 
GP(m_^i) = amplitude of the concentrate vo^tex equivalent to 
GB(m^i) . 

GL(m^i) = amplitude of the concentrate vortex equivalent to 
GR(m^i). 

The integral appearing in eq. (43) is evaluated analytically 



^m+1 

J — = im. I2(D+ p- r cos(cp(i) + - 0)) 






The value of the bound circulation is obtained either directly 
if the loading is expressed as in eq. (21) or after solving a Volterra 
equation (25) if it is expressed as in eq. (20). 

The solution of the Volterra equation can be obtained using the 
Liouvil le-Neumann method of successive substitutions [20]. This gives 



= o'-'H'rT + ■ 

^ t o 

CO (ik -ik(r) 

Tj (-1) ' cos 0 e • \ Ap(r , t) sintdt] 

k=1 '^0 



( 46 ) 
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When the pressure loading is given in the form of Birnbaum modes 
(20)^ the above integral is evaluated using the relations 



\ t /V /V 

\ cotan sin t dt = 0 - sin0 




/\ 

\ sin nt sin t dt = it - ^ sin(20) , 


n = 1 


sin((n - 1 ) *0) 
2(n-1) 


sin((n +1 ) -0) , 

2(n+1) ^ n 1 



i 



{ 



i 



( 

I 
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Chapter IV 
RESULTS 

The numerical solution developed in the preceding chapter was 
put into computer form; the details of the program are presented in 
Appendix B. In order to test the workability of this program we 
used it to compute the pressure at an upstream field point using the 
NSRDC Propeller No. 4ll8 [22], The design characteristics and oper-* 
ating conditions of this propeller are represented in Figs. 7 and 8. 
We considered both the open water condition and the nonuniform wake 
obtained with a screen which produces a circumferentially varying 
sinusoidal wake composed only of multiplies of blade harmonics. 
Unfortunately^ we did not have information about the loading distri- 
bution for those conditions and therefore we used an arbitrary number 
of modes of the Birnbaum distribution with coefficients equal to 
unit. In order to test the solution of the Volterra equation (46) 

the input loading was expressed in terms of pressure distribution 

\ 

rather than vorticity distribution. 

We made several sample calculations in order to evaluate the 
effect of changing the grid spacings and to determine how far down- 
stream we needed to carry the streamwise integrations in the wake 
(Table 1 shows a sample of the computer output for one of the cal- 
culations). We found that for the field point considered (x = “*8^ 
r = 1 . 2^ 6 = 0) ^ we needed about 20 chordwise panels and had to keep 
the angular grid spacing down to 2 degrees. The contribution due to 



m 

\ 
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r/R— 1 • 0 




r/R 


1/D 


T/1 


P/D 


0.2 


0.3200 


0.1028 


1.0860 


0.3 


0.3635 


0.0776 


1.0845 


0.4 


0.4048 


0.0590 


1.0823 


0.5 


0.4392 


0.0451 


1.0796 


0.6 


0.4610 


0.0348 


1.0780 


0.7 


0.4622 


0.0271 


1.0766 


0.8 


0.4347 


0.0210 


1.0750 


0.9 


0.3613 


0.0166 


1.0735 


1.0 






1.0700 



FIG. 7 - PROPELLER 4118 DESIGN CHARACTERISTICS 
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FIG. 8 - RADIAL DISTRIBUTION OF THE AMPLITUDES AIJD 
PHASES OF THE THIRD HARKONIC OF THE 
CIRCUNFERENTIAL WAKE VARIATION 
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PRESSURE FIELD AROUND A MARINE PROPELLER 



PROPELLER CHARACTERISTICS 

TYPE NSRDC, PROPELLER NO . 4118 

DIAMETER 12.0000 INCHES 

HUB DIAMETER 2.4000 INCHES 

NO. OF BLADES 3 

RPM 1080.0000 



BLADE OUTLINE (UNITS -FRACTION OF RADIUS) 



RADII 


LEAD. EDGE 


TRAIL. EDGE 


MAX.TH 


0.250 


-0.3417 


0.3417 


0.0605 


0.300 


-0.3635 


0.3635 


0.0564 


0.400 


- 0.4048 


0.4048 


0.0478 


0.500 


-0.4392 


0.4392 


0.0396 


0.600 


- 0.4610 


0.4610 


0.0321 


0.700 


- 0.4615 


0,4615 


0.0250 


0.800 


-0.4347 


0.4347 


0.0183 


0.900 


-0.3613 


0.3613 


0.0120 


0.950 


-0.2500 


0.2500 


0.0065 


1.000 


0.0 


0.0 


0.0 



GRID SPACING 

NO. OP CHORDWISE PANELS 20 

ANGULAR CHORDWISE SPACING 0.0349066 RADIANS 



TABLE 1 - SAMPLE COMPUTER OUTPUT 
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inflow CHARACTERISTICS 



SYMBOLS AND UNITS 

VS SHIP SPEED, KNOTS 

RADII FRACTION OF RADIUS 

XVXO ZEROTH HARMONIC OF INFLOW / VS 

TANBETA TANGENT OF ADVANCE ANGLE 

NHARM HARMONIC NUMBER 

XVXN AI4PLITUDE OF NHARM OF INFLOW / VS 
XPHN PHASE OF NHARM OF INFLOW, RADIANS 



vs= 


8.8780 


NHARH= 3 






RADII 


XVXO 


TANBETA 


XVXN 


XPHN 


0.250 


1.0000 


1.0606 


0.0800 


- 0.2269 


0.300 


1.0000 


0.8838 


0.1000 


- 0.1920 


0.400 


1.0000 


0.6629 


0.1500 


-0.1047 


0.500 


1.0000 


0.5303 


0.2000 


0.0 


0.600 


1.0000 


0 . 4419 


0.2500 


0.0873 


0.700 


1.0000 


0.3788 


0.2900 


0.0 


0.800 


1.0000 


0.3314 


0.3200 


-0.0314 


0.900 


1.0000 


0.2946 


0.3400 


0.0 


0.950 


1.0000 


0.2791 


0.3450 


0.0436 



TABLE 1 (cent.) 
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loading DISTRIBUTION 
(VORTICITY MODES) 



RADII 


MODE 


REAL 


0.250 


1 


0.014860 


0.250 


2 


0.014860 


0.250 


3 


0.014860 


0.300 


1 


0.014860 


0.300 


2 


0.014860 


0.300 


3 


0.014860 


0,400 


1 


0.014860 


0.400 


2 


0.014860 


0.400 


3 


0.014860 


0.500 


1 


0.014860 


0.500 


2 


0.014860 


0.500 


3 


0.014860 


0.600 


1 


0.014860 


0.600 


2 


0,014860 


0.600 


3 


0.014860 


0.700 


1 


0.014860 


0.700 


2 


0.014860 


0.700 


3 


0.014860 


0.800 


1 


0.014860 


0.800 


2 


0.014860 


0.800 


3 


0.014860 


0.900 


1 


0.014860 


0.900 


2 


0.014860 


0.900 


3 


0.014860 


0.950 


1 


0.014860 


0.950 


2 


0.014860 


0.950 


3 


0.014860 



COEF. 



IMAG. 

0.014860 

0,014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 

0.014860 



TABLE 1 (CONT) 
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PRESSURE FIELD RESULTS 



SYMBOLS AND UNITS 

EX AXIAL COORDINATE OF FIELD POINT, 

FRACTION OF RADIUS 

ER RADIAL COORDINATE OF FIELD POINT, 

FRACTION OF RADIUS 

TETA ANGULAR COORDINATE OF FIELD POINT, 
RADIANS 

NT NONDIMENSIONAL TIME, NT=TIME*RPM/60 

KP PRESSURE COEFF., P/ ( DENSITY* ( (RPM/ 

60)**2)*(DIAMETER**2) ) 



EX= -0.8000 ER= 1.2000 TETA= 0.0 

1 - TOTAL PRESSURE 

NT KP 

REAL IMAG 



0.0 


0.109260 


0.064822 


0.0625 


-0.002797 


O.O87227 


0.1250 


- 0.101459 


0.032849 


0.1875 


-0.086291 


-0.053374 


0.2500 


0.115199 


-0.156751 


0.3125 


O.O85257 


-0.032765 


0.3750 


0.021487 


0.064308 


0.4375 


-0.051879 


0.049902 


0.5000 


-0.0304070 


- 0.024493 


0.5625 


-0.045830 


- 0.065644 


0.6250 


0.036785 


0.004191 


0.6875 


0.028781 


0.083445 


0.7500 


-0.047159 


0.067708 


O.8I25 


- 0.08040 


-0.061579 


0.8750 


0.040268 


-0.055254 


0.9375 


0.117604 


0.020541 



TABLE 1 (cont) 
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2 - HARMONIC COMPONENTS 
HARM. NO. KP 





REAL 


IMAG 


0 


0.000563 


0.000171 


1 


-0.000115 


0.001645 


2 


0.008606 


0.000716 


3 


0.069201 


0.047710 


4 


-0.001695 


-0.002692 


5 


-0.011514 


-0,011500 


6 


-0.010340 


0.001098 


7 


0,000716 


0.001721 


8 


0.010076 


0.001721 


9 


-0.001120 


-0.009159 


10 


-0.000015 


0.000259 


11 


0.000433 


0.000063 


12 


0.019964 


-0.009159 


13 


0.000831 


-0.000075 


14 


-0.000070 


-0.000130 


15 


0.018597 


0.007825 



TABLE 1 (cont.) 
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the wake after about 7 revolutions was found negligible^ which agrees 
with calculation by Kerwin [18] for the zeroth harmonic* For field 
points down stream the wake integration is carried down 7 revolutions 
beyond the axial position of the field point. 

It is shown in Table 2 how the number of chordwise panels and 
the angular grid spacing affect the value of the harmonic components 
of the pressure. 

Figures 9 and 10 show the pressure variation and harmonic 
components at the field point due to the zeroth and third harmonic 
of loading respectively. 

Since we used an arbitrary loading' distribution we cannot assess 
the accuracy of the results. This will have to be done at a later 
date. 
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B 



No, of panels 
Grid spacing 



15 

4 degrees 



15 

2 degrees 



HARM. NO. Kp 

REAL ^ IMAG 

0 0.0 0.0 

3 0.081273 0.054870 

6 -0.014081 0.000345 

9 0.002288 -0.017413 



K, 



REAL 

0.0 

0.076421 



IMAG 

0.0 

0.053042 



•O.OO9OI8 0.000736 

0.001631 -0.012273 



No. of panels 
Grid spacing 



20 

2 degrees 



HARM .NO. 

0 

3 

6 

9 



K, 



REAL 

0.0 



-IMAG 

0.0 



0.069201 0.047710 
•0.010340 0.001098 
0.001120 -0.009159 



TABLE 2 - EFFECT OF NUMBER OF CHORDWISE PANELS 



AND AIJGULAR GRID SPACING 
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fig. 9. a - PRESSURE AT UPSTREAM FIELD POINT DUS TO ZEROTH HARHIONIC OP LOADING 
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Field point: x/R = 0.8 

r/R = 1.2 

0 = 0 , 

HARMONIC COMPONENTS 
HAPiM.NO, KP 





real 


imag. 


0 


-0.059293 


0.0 


1 


- 0.000003 


-0.000004 


2 


-0.000002 


-0.000002 


3 


-0.001161 


-0.064603 


4 


-0.000003 


-0.000001 


5 


0.000000 


-0.000001 


6 


-0.000897 


-0.000608 


7 


-0.000006 


0.000004 


8 


0.000000 


0.000000 


9 


-0.000006 


-0.000004 



FIG. 9.b - PRESSURE AT UPSTREM FIELD POINT DUE TO ZEROTH 



HARMONIC OF LOADING - HARI-IONIC COMPONENTS 
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I 



h 



real part 
imaginary part 
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10, a - PRESSURE at upstream FIELD POINT DUS TO THIRD HARI.10NIC OF LOADING 
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Field point: x/R = 0.8 

r/R = 1.2 

0 = 0 . 

.HARMONIC COMPONENTS 
HARM. NO. KP 



0 


real 

0^000563 


lEiag. 

o;oooi 7 i 


1 


-0.000115 


0.001645 


2 


0.008606 


0.000716 


3 


0.069201 


0.047710 


4 


-0.001695 


-0.002692 


5 


-0.001514 


-0.001500 


6 


-0.010340 


0.001098 


7 


0.000716 


0.001721 


8 


0.000076 


-0.000898 


9 


- 0.001120 


-0.009159 



FIG. 10. b - PRESSURE AT UPSTREAM FIELD POINT DUE TO 
THIRD HARMONIC OF LOADING - HARMONIC 



COMPONENTS 
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Chapter V 

CONCLUSIONS AND RECOMMENDATIONS 

A computerized procedure has been developed for the calculation 
of the pressure field generated by a propeller operating in an un- 
steady flow. We believe that it represents an important improvement 
over the existing methods based on lifting line theory and circum- 
ferentially uniform wake. It is also a useful complement to the 
recently published work by Jacobs ^ [19] which develops a 

different approach to the same problem. 

A considerable amount of time was put into the optimization of 
the program but it is probable that it can still be ameliorated. 

One of the most time consuming operations -- the integration of the 
streamwise vorticity in the wake -- was successfully reduced by the 
use of the Euler Transformation. The two other operations whose 
improvement would have an important effect on the over -all improvement 
of the program are the integration of the radial vorticity in the 
wake and the calculation of the velocity potential. 

Sample calculations were made to determine the sensibility of 
the results to the number of chordwise panels^ angular grid spacing 
and how far downstream integrations were carried. Assessment of the 
accuracy of the results was not possible since v/e used an arbitrary 
loading distribution. 

Certainly^ the next thing that needs to be done is to compare 
results with experimental data or with the results obtained in [19], 






This requires that the loading distribution be obtained by any 
of the existing methods referred to before. The program accepts 
the input loading distribution in terms of pressure or circulation. 
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Appendix A 



VELOCITY POTENTIAL INDUCED BY A SINGLE VORTEX 



In this appendix we compare the velocity potential induced by 
a single reentrant vortex with the one induced by a concentrated 
doublet of equivalent strength. 

The velocity potential due to a vortex is given by 



S = vortex strength 
A - area enclosed by the vortex 

Y = angle between the normal to A and the line from A 
to the field point 

D = distance from A to the field point. 



Equation (A-1) is similar to the one that gives the velocity 
potential induced by a uniform distribution over A of source doublets 
of strength S. It is clear that the strength of the equivalent con- 
centrated doublet must be equal to S x A, 

In order to simplify the calculations^ we will consider a 
rectangular vortex centered at the origin^ as shown in Fig, A-1. 

The velocity potential at p is found using eq, (A-1) 




(A-1) 



where 
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FIG. A-1 - ORIENTATION OF REENTRANT VORTEX 



The above integral can be evaluated analytically by integrating 
directly in one of the variables -- say -- and then making the 
following changes of variables before performing the next integration: 

71 - y 

[(x a)^ + 

2 - 

u=v+(v +1)^ 

2 

t = u 

The final result is 

$v == [(tan^A - tan^B) - (tah^C - tan“D)] (A-2) 

where 

A = (y+ a)^ + z^ - (y - a)[(x+ a)^+ (y - a)^+ 

z(x + a) 
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B = 



c = 



D = 



(y+a)^+z^-(y + a)f(x + a)^ + (y + a)^+z^]^ 
z(x +a) 

(y - a)^ +z^ - (y - a)r (x - a)^ + (y - a)^ + z^]^ 
z(x - a) 

(y + a)^ + z^ - (y + a)[(x~a)^+(y + a)^ + z^]^ 
z (x - a) 



The velocity potential induced at p by a concentrated doublet 



of strength 4a located at the origin is given by 



^ a^Sz 

^d ” , 2 . 2 . 2,^ 



(A- 3) 



7T(x +y +z 

It is physically reasonable to predict that the largest dif- 
ferences between the values given by eqs. (A-1) and (A-2) will be found 
at the field points lying on the z axis. Figure A-2 shows how those 
differences vary with the distance along the z axis. 

It is apparent from the analysis of Fig. A-2 that the ratio of 
the two values approaches very rapidly the unity. For z/a = 8 (i.e.^ 
for a distance 4 times greater than the side of the original vortex) 
the difference between the two values is less than one percent. 
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Appendix B 



COMPUTER PROGRAM 



1 . General 

The program presented in this appendix is essentially the 
translation into FORTRAN IV of the numerical solution outlined in 
Chapter III. 

The initial section of the program in which the blade and wake 
are divided into panels is taken from Professor Kerwin’s program^ [18], 
Also^ the scheme used for the zeroth harmonic of the loading is based 
on Professor Kerwin*s work. 

The program is written for use on the IBM Systenv/360 digital 
computer. For different processors slight changes might have to be 
made. The amount of memory space required is about 150K bytes and 
the computation time for a single field pointy including harmonic 

i 

analysis of the pressure is about % minutes for the zeroth harmonic 
and S minutes for higher harmonics^ using a FORTRAN G Compiler. 

The complex harmonic analysis is performed by the subroutine 
HARM described in [23]. The final output variables are subscripted. 
This was done in anticipation of a plotting subroutine which can be 
added in the future. 
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2 . List of Key Variables 

The following is a list of the principal variables in the program^ 



arranged 


in alphabetical order. 


A(M,N) 


= Coefficients of the Birnbaum modes of the loading function 


BOUNDA 


= Axial velocity of the field pointy induced by the bound 
vortices 


BOUND! 


= Tangential velocity at the field pointy induced by the 
bound vortices 


BWAKEA 


= Axial velocity at the field pointy induced by the radial 
vortices in the wake 


BWAKET 


= Tangentical velocity at the field pointy induced by the 
radial vortices in the wake 


CHOR(M) 


= Local chord at radius R(M) 


DHUB 


= Hub diameter 


DIAM 


= Propeller diameter 


ER 


= Radial position of the field point 


EX 


= Axial position of the field point 


G(M) 


= Chordwise integral of the bound circulation at radius R(M) 


GB(M,N) 


= Bound circulation assumed concentrated on the N^th grid 
line in the M^th chordwise panel 


GL(M,N) 


= Wake circulation assumed concentrated between the 

KWAKE(M) + N and KWAKE(M) + N + 1 grid lines in the M'th chord 
wise panel^ for computing PTWK 


GLR(M) 


= Advance coefficient at radius R(M) 



o 






GLS 

GP 

GR 

GT(M,N) 

IDENT 

IGB 



KD(N) 

KN(M) 

KWAKE(M) 

MN 

NBLADE 

NHARM 

NL 

NLE(M) 

NMAX 

NR 

NT 

NTE(M) 



Advance coefficient based on VS 

Circulation in the blade^ assumed concentrated in the 
middle of each element of the lattice for computing PTBL 
Radial circulation in the wake 

Streamwise circulation in the blade^ between the N and 
N + 1 grid lines in the M’th panel 
Identification of the propeller 

Parameter that describes the loading function. If IGB = 0 
the loading is in terms of bound circulation; if IGB > 0 
it is in terms of the pressure jump in the blade 
Angle between the N’th blade and the YY axis 
Reduced frequency at radius R(M) 

Identification of the first grid line of the wake ar 
radius RZ(M) 

Number of equal spaces in the division of the span into 
chordwise panels 
Number of blades 
Harmonic number 

Largest integer multiple of the grid spacing with value 
less or equal to TMIN 

Identification of the first grid line in the M^th panel 
Number of loading modes 

Largest integer multiple of the grid spacing with value 
less or equal to TMAX 

Total number of grid lines in the blade 

o 

Identification of the last grid in the blade in the M'th panel 
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NTH 


= Parameter that defines the spacing between grid lines as 
NTH X 2 degrees 


NW(M) 


= Total number of grid lines spanned by the blade in the 
M*th panel 


PHN (M) 


= Phase angle of the NHARM*th harmonic of the inflow velocity 


PRS (M) 


= Field point pressure when the first blade is at an angle 
of M*7 i/( 8 NBLADE) with fixed vertical axis 


PTBL 


= Time derivative of the velocity potential at the field 
pointy induced by the circulation in the blade 


PTWK 


= Time derivative of the velocity potential at the field 
pointy induced by the circulation in the wake 


PTSC 


= Time derivative of the velocity potential at the field 
pointy induced by the source distribution 


R(M) 


= Midradius of the M* th chordwise panel 


RPM 


= Rotation per minute 


RZ(M) 


= Radius at the lower bound of the M*th panel 


SB(M,N) 


= Source strength concentrated on the N*th grid line in 
the M*th chordwise panel 


SL(M) 


= Leading edge at radius R(M) 


SMA(N) 


= Differential of the thickness at the N*th grid line^ as 
percent of the max thickness in each panel 


SOURCA 


= Axial velocity at the field pointy induced by the source 
distribution 


SOURCT 


= Tangential velocity at the field pointy induced by the 
source distribution 
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ST(M) 


= Trailing edge at radius R(M) 


STAR(M) 


= Factor which multiplied by SMA(N) gives the source strength 
at the N^th grid line in the M*th chordwise panel 


TETA 


= Angular coordinate of the field point 


TK(N) 


= Standard thickness^ fraction of the max thickness 


TL(M) 


= Amplitude of the streamwise vorticity in the wake at 
radius RZ(M) 


TLE (M) 


= Angular coordinate of the leading edge at radius R(M) 


TMAX 


= Extreme value of TTE(M) over all the panels 


TMIN 


= Extreme value of TLE(M) over all the panels 


TRAILA 


= Axial velocity at the field pointy induced by the 
streamwise vorticity in the blade 


TRAILT 


= Tangential velocity at the field pointy induced by the 
streamwise vorticity in the blade 


TTE (M) 


= Angular coordinate of the trailing edge at radius R(M) 


TWAKEA 


= Axial velocity at the field pointy induced by the stream- 
wise vorticity in the wake 


TWAKET 


= Tangential velocity at the field pointy induced by the 
streamwise vorticity in the wake 


VFP 


= Inflow velocity at the field point 


VS 


= Ship speed 


VXO(M) 


= Zeroth harmonic of the inflow velocity at radius R(M) 


XGL(N) 


= Advance coefficient at radius XR(N) 


XKN (N) 


= Reduced frequency at radius XR(N) 
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XPHN(N) 

XR(N) 

XSL(N) 

XST(N) 

XSTAR(N) 

XTBETA(N) 

XTZ(N) 

XVXN(N) 



= Phase angle of the NHARM' th harmonic of the inflow 
velocity at radius XR(N) 

= Input radius (see Input Description) 

= Leading edge at radius XR(N) 

= Trailing edge at radius XR(N) 

= Factor of proportionality of the source strength at 
radius XR(N) 

= Tangent of the advance angle at radius XR(N) 

= Maximum thickness at radius XR(N) 

= NHARM' th harmonic of the inflow velocity at radius XR(N) 
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3 . Input Description 

The following sequence of input cards is required for the 
operation of the program: 

A IDENTIFICATION 1 card 

B GENERAL DESCRIPTION 1 card 

C BLADE OUTLINE 3 cards 

D WAKE CHARACTERISTICS 4 cards 

E LOADING DISTRIBUTION 21 cards max. 

F FIELD POINT DESCRIPTION 1 card 

A complete description of the above input cards is presented in 
Table B-1. The input radii referred to -in that table are schematically 
shown in the figure below. 

XR(1) XR(2) XR(9) XR(IO) 

Yi I I II \\\^ 

r^/R 1.0 

DX 

FIG. B-1 - INPUT RADII OF PROPELLER DATA 



DX is the Input Spacing^ equal to ^ X interval from hub to tip. 
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